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ABSTRACT 

A  queuing  system  is  investigated  in  which  two  types 
of  customers,  Type  1  and  Type  2,  attempt  to  access  a  multi- 
server  operator  switchboard.   Customer  arrivals  for  each 
type  are  assumed  to  be  in  accordance  with  independent 
Poisson  processes  with  time  constant  parameters.   Service 
times  are  assumed  to  be  exponentially  distributed  and  each 
type  of  customer  has  a  different  service  rate.   When  a 
queue  is  formed,  Type  1  customers  are  allowed  to  queue 

regardless  of  the  length  of  the  queue  but  Type  2  customers 
are  blocked  when  the  queue  size  exceeds  a  predesignated 
limit. 

An  exact  solution  of  the  problem  was  not  obtained  and 
a  discussion  of  the  complications  which  precluded  an  exact 
solution  is  presented.   A  simulation  algorithm  is  developed 
and  validated  by  comparison  with  an  exact  solution  of  the 
special  case  of  the  problem  in  which  the  service  rates 
are  equal. 

Finally,  an  analytical  approximation  is  developed  and 
shown  to  give  reasonable  results  by  comparison  with  the 
simulation  model. 
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I.   INTRODUCTION 

A.   STATEMENT  OF  THE  PROBLEM 

The  problem  addressed  in  this  thesis  arises  in  tele- 
communications systems  and  concerns  the  flow  of,  and  delays 
to,  customer  traffic  through  a  multiple  operator-serviced 
switchboard.   There  are  two  types  of  customers,  Type  1  and 
Type  2,  and  each  type  has  its  own  arrival  and  service  parameters 
The  arrival  and  service  processes  of  the  two  types  are  sta- 
tistically independent.   All  switchboard  operators  can  ser- 
vice either  type  of  customer,  and  the  queuing  discipline  is 
first- in,  first-out  (FIFO) .   What  makes  this  problem  inter- 
esting is  that  while  Type  1  customers  are  allowed  to  queue 
regardless  of  the  queue  size,  Type  2  customers  are  blocked 
from  joining  the  queue  when  the  queue  length  equals  or  exceeds 
a  predesignated  limit. 

An  example  might  be  useful  to  ensure  that  the  problem  is 
understood.   Consider  the  civil-sector  telephone  system,  and 
let  the  Type  1  customers  be  those  who  are  direct-dialing  a 
collect,  long-distance  call,  and  the  Type  2  customers  be 
those  who  wish  the  operator  to  dial  the  call  for  them.   The 
two  types  can  be  differentiated  by  the  system  on  the  basis 
of  the  number  of  digits  they  dial;  the  Type  1  customers  would 
dial  a  total  of  eleven  digits  while  the  Type  2  customers 
would  dial  the  single  digit  0.   Since  the  Type  1  customers 
require  less  work  on  the  part  of  the  operator,  it  may  be 


desirable  under  conditions  of  heavy  use  when  the  queue  length 
gets  large,  to  discourage  customers  from  having  the  operator 
dial  their  calls  for  them  by  blocking  Type  2  customers  from 
joining  the  queue.   In  this  case,  both  types  of  customers 
would  be  allowed  to  queue  when  the  queue  length  is  small, 
perhaps  receiving  a  recording  advising  them  that  the  operators 
are  busy  and  asking  them  to  wait.   However,  when  the  queue 
length  exceeds  its  predesignated  limit,  Type  2  customers 
receive  a  busy  signal  or  possibly  a  recording  advising  them 
to  direct-dial  their  call. 

1.   Simplifications  and  Assumptions 

In  order  to  make  the  problem  tractable,  a  number  of 
complications  which  would  exist  in  a  real-world  system  have 
been  ignored  and  a  number  of  simplifying  assumptions  have 
been  made.   Among  those  elements  of  the  problem  which  are 
not  being  considered  are  the  following: 

a.  Preemption  by  a  high  priority  class  of  customers. 
A  particular  example  would  be  customers  who  dial  a  three- 
digit  code  in  the  event  of  an  emergency  and  who  must  receive 
preferential  service. 

b.  Balking  or  reneging.   For  purposes  of  this 
study,  it  is  assumed  that  once  a  customer  has  joined  the 
queue,  he  remains  there  until  his  service  commences. 

c.  The  existence  of  a  retrial  population.   In  an 
actual  situation,  a  Type  2  customer  who  is  blocked  from  queuing 
would  be  expected  to  do  one  of  three  things:   (i)  decide  not 

to  attempt  to  complete  the  call  at  this  time;  (ii)  try  accessing 


the  system  as  a  Type  1  customer;  or  (iii)  continue  trying 
to  access  the  system  as  a  Type  2  customer  at  repeated,  random 
intervals.   The  repeated  attempts  to  access  the  system  by 
those  Type  2  customers  who  exercise  option  (iii)  can  be  seen 
to  have  an  effect  on  the  arrival  rate,  and  these  customers  are 
said  to  belong  to  the  retrial  population.   Although  there 
exist  a  number  of  ways  to  model  this  phenomena,  no  attempt 
has  been  made  to  do  so  in  this  study.   In  addition,  the 
increase  in  the  arrival  rate  of  Type  1  customers  caused  by 
those  Type  2  customers  who  choose  option  (ii)  is  ignored. 
The  assumptions  which  were  considered  necessary  were: 

d.  Arrivals  are  in  accordance  with  a  Poisson  process 
and  the  arrival  rates  are  constant  over  time  for  each  type 

of  customer.   This  is,  admittedly,  an  unrealistic  assumption 
since  it  is  unlikely  that  the  customer  arrival  rate  will  be 
the  same  at  two  A.M.  as  it  is  at  two  P.M.   Although  arrivals 
according  to  an  inhomogeneous  Poisson  process  can  be  modeled 
fairly  well  by  computer  simulation,  the  introduction  of  time- 
dependent  arrival  rates  greatly  decreases  the  chances  of 
finding  an  analytical  solution  or  approximation  and  it  is 
for  this  reason  that  they  have  not  been  considered. 

e.  Service  times  are  exponentially  distributed  and 
service  rates  are  constant.   This  assumption  is  somewhat  less 
unrealistic  since  it  seems  plausible  that  service  times  should 
be  the  same,  on  the  average,  regardless  of  the  time  of  day. 


2.   Implications  of  Unequal  Service  Rates  and  the 
Queuing  Discipline 

The  facts  that  the  two  types  of  customers  have 

different  service  rates  and  that  Type  2  customers  are  blocked 

from  the  queue  when  it  becomes  too  long  can  be  seen  to  have 

several  important  implications  for  the  way  in  which  the  system 

functions.   In  particular,  consider  the  following: 

a.  Since  the  service  rates  are  unequal,  the  system 
service  rate  is  dependent  not  only  on  the  number  of  custo- 
mers in  service,  but  also  on  their  types.   Therefore,  any 
solution  technique,  either  analytical  or  simulation,  must  be 
able  to  keep  track  of  the  composition,  by  type,  of  the  group 
of  customers  undergoing  service  at  every  time. 

b.  When  a  change  occurs  in  the  state  of  the  system, 
the  new  composition  of  the  group  of  customers  undergoing 
service  depends  not  only  on  whether  there  are  customers       -. 
waiting  in  the  queue,  but  also  on  the  precise  type  of  the 
customer  at  the  head  of  the  queue.   Therefore,  it  is  necessary 
to  keep  track  of  the  composition,  by  type,  of  the  queue,  in 
detail. 

c.  Since  the  blocking  of  Type  2  customers  effectively 
sets  their  arrival  rate  to  zero  when  the  condition  occurs, 

it  can  be  seen  that  the  overall  system  arrival  rate  is  not 
constant,  but  is  dependent  upon  the  number  in  the  system. 
This  effect  produces  an  interesting  problem  in  computing 
the  expected  waiting  time  for  an  arriving  customer.   This 
problem  is  discussed  in  Chapter  Three. 


10 


3 .   Figures  of  Merit 

Three  figures  of  merit  were  chosen  to  serve  as  indi- 
cators of  system  performance.   They  are:   (1)  the  expected 
waiting  time  for  a  customer  of  either  type  who  joins  the 
system,  (2)  the  proportion  of  time  that  there  are  no  customers 
of  either  type  in  the  system  and  (3)  the  proportion  of  time 
that  the  queue  length  equals  or  exceeds  its  limit  and  Type  2 
customers  are  blocked.   For  all  three  of  the  above,  only  the 
steady  state  or  long-run  value  is  considered. 

B.   ANALYTICAL  PROCEDURE 

1.   The  Birth-Death  Process 

The  first  attempt  at  a  solution  was  to  model  the 
system  as  a  multi-dimensional,  birth-death  process  and  then 
attempt  to  solve  the  resulting  balance  equations  algebraically. 
The  necessity  of  keeping  track  of  the  composition,  by  type, 
of  both  the  group  of  customers  undergoing  service  and  the 
group  of  customers  in  the  queue  meant  that  the  minimum  number 
of  dimensions  which  could  be  used  to  adequately  model  the 
state  of  the  system  was  four.   That  is,  the  state  variable  which 
represented  the  total  number  of  customers  in  the  system  would 
have  four  dimensions;  one  each  for,  the  number  of  Type  1 
customers  undergoing  service,  the  number  of  Type  2  customers 
undergoing  service,  the  number  of  Type  1  customers  in  the 
queue  and  the  number  of  Type  2  customers  in  the  queue.   How- 
ever, even  using  a  four  dimensional  state  variable  proved  to 
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be  inadequate  to  differentiate  between  the  system  state  in 
which  a  Type  1  customer  was  first  in  the  queue  and  the  state 
in  which  a  Type  2  customer  was  the  first  in  the  queue.   Since 
the  way  in  which  the  system  changes  states  is  dependent  upon 
the  type  of  the  customer  at  the  head  of  the  queue,  it  would 
have  been  necessary  to  expand  the  state  variable  beyond  four 
dimensions  to  model  this  effect.   The  complexity  of  the 
balance  equations  in  the  four  dimensional  case  offered  little 
hope  of  finding  a  solution  and  an  expansion  beyond  four  dimen- 
sions would  reduce  the  likelihood  of  solving  the  equations 
still  further.   For  this  reason,  the  multi-dimensional,  birth- 
death  process  approach  was  abandoned  and  other  alternatives 
were  considered. 

2.  Simulation 

The  next  attempt  to  model  the  system  was  via  computer 
simulation  since  it  was  hoped  that  the  development  and  exe- 
cution of  the  simulation  model  would  provide  insights  into 
the  workings  of  the  system  which  would  lead  to  ideas  for 
analytical  solutions  or  approximations.   The  development  of 
the  simulation  model  is  covered  in  detail  in  Chapter  two. 

3.  Verification  of  the  Simulation  Model 

Although  verification  of  a  simulation  model  is  a 
normal  step  in  the  development  of  the  model,  it  is  mentioned 
separately  here  for  several  reasons.   First,  it  was  a  major 
effort  in  itself;  secondly,  it  utilized  a  generalization  of 
Little's  equation  for  the  waiting  time  that  might  be  of  some 
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interest;  and  finally,  it  served  as  the  basis  for  the  first 
attempt  at  an  analytical  approximation.   The  verification 
of  the  simulation  is  covered  in  detail  in  Chapter  Three. 
4 .   An  Analytical  Approximation 

Chapter  Four  details  the  derivation  of  an  approximate 
analytical  solution  based  on  the  equations  developed  in 
Chapter  Three  for  verifying  the  simulation  model.   In  addition, 
the  rationale  for  developing  an  analytical  approximation 
model  and  the  advantages  to  be  gained  from  using  such  a  model 
are  presented.   Finally,  comparisons  are  made  between  results 
obtained  from  the  analytical  approximation  and  the  simulation 
to  determine  the  degree  of  agreement  between  the  two  models . 


13 


II.   THE  SIMULATION  ALGORITHM 

A.   RATIONALE  FOR  THE  SIMULATION  APPROACH 

Although  it  is  generally  accepted  that  an  analytical 
solution  is  preferable  to  a  simulation,  there  are  a  number 
of  situations  in  which  a  simulation  approach  is  indicated. 
Foremost  among  these  is  the  case  in  which  an  analytical 
solution  does  not  exist  or,  equivalently ,  an  analytical  solu- 
tion exists  in  theory,  but  is  unobtainable  in  practice.   As 
was  mentioned  in  Chapter  One,  the  latter  situation  is  the  case 
in  the  present  problem.   If  it  is  impractical  to  obtain  a 
complete  and  useable  analytical  solution,  a  search  for  a  simple 
analytical  approximation  is  likely  to  be  worthwhile.   Given 
the  complexities  of  the  problem  of  this  thesis,  derivation  of 
an  analytical  approximation  becomes  a  non-trivial  task. 
Therefore,  the  simulation  approach  is  utilized  to  gain  insight 
into  the  problem  and  the  operation  of  the  queuing  system  in 
order  to  facilitate  later  analytical  work. 

Another  reason  for  simulating  the  system  is  that  if  an 
analytical  approximation  can  be  found,  then  the  simulation 
can  serve  as  a  basis  of  comparison.   This,  of  course,  assumes 
a  high  degree  of  confidence  in  the  simulation  model  which 
means  that  the  simulation  must  somehow  be  verified.   This 
important  step  in  the  development  of  the  simulation  model  is 
covered  in  Chapter  Three. 
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Finally,  if  a  suitable  analytical  approximation  cannot 
be  found,  then  the  simulation  remains  the  only  tool  with 
which  to  study  the  system. 

B.   CONCEPTUAL  OPERATION  OF  THE  SYSTEM 

The  first  step  in  the  development  of  the  simulation 
algorithm  was  to  study  the  operation  of  the  system  so  as 
to  make  sure  that  the  manner  in  which  the  system  changed 
from  one  state  to  another  was  fully  understood.   The  system 
was  viewed  as  staying  in  a  current  state  for  a  random  time 
after  which  one  of  four  statistically  independent  events  would 
occur  with  certain  specified  probabilities.   The  occurrence 
of  these  events  causes  the  system  to  change  state.   The  four 
possible  events  are:   (i)  a  Type  1  customer  arrives  and  joins 
the  system,  (ii)  a  Type  2  customer  arrives  and  may  or  may 
not  join  the  system  depending  on  the  current  state,  (iii)  a 
Type  1  customer  completes  service,  or  (iv)  a  Type  2  customer 
completes  service.   The  conceptual  operation  of  the  system 
and  the  four  events  along  with  their  probabilities  of 
occurrence  are  displayed  in  Figure  2.1.   Figure  2.2  shows  the 
system  changes  associated  with  the  occurrence  of  each  of  the 
events.   Terms  used  in  the  two  figures  are  defined  below. 

Jl   =   number  of  Type  1  customers  in  service. 

J2   =   number  of  Type  2  customers  in  service. 

Ql   =   number  of  Type  1  customers  in  the  queue. 

Q2   =   number  of  Type  2  customers  in  the  queue. 

S    =   number  of  servers. 
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T   =   queue  threshold  (applies  to  Type  2  customers) . 

N    =   total  number  of  customers  in  the  system 
=  Jl  +  J2  +  Ql  +  Q2 

X        =   Type  1  arrival  rate. 

A 2   =   Type  2  arrival  rate. 

y,   =   Type  1  service  rate. 

y2   =   Type  2  service  rate. 

d    =  X1    +    X2    +   Jly   +  J2y2 
A  few  notes  will  be  helpful  to  fully  understand  the  two 
figures. 

1.   The  "sojourn  time"  is  the  random  time  period 
during  which  the  system  stays  in  its  current  state.   The 

sojourn  time  is  exponentially  distributed  with  mean  equal  to 

i  .       i 

which  is  -=■•   Note  that  this  means  that 


X      +    X2    +   Jly   +  J2y2      w"     d 

the  distribution  of  the  sojourn  time  is  not  fixed,  but  that 
it  changes  with  each  change  in  Jl  and/or  J2 . 

2.  According  to  the  queuing  discipline  assumed  for 
this  problem,  a  Type  2  customer  can  only  join  the  system 
when  N,  the  number  in  the  system,  is  less  than  T,  the  queue 
threshold.   This  is  illustrated  in  Figure  2.2,  where  it  can 
be  seen  that  no  state  change  takes  place  if  a  Type  2  arrival 
occurs  when  the  number  in  the  system  is  larger  than  the  cutoff 
value. 

3.  Figure  2.2  also  illustrates  the  fact  that  when 
a  queue  is  formed,  the  new  state  of  the  system  depends  upon 
whether  a  Type  1  or  Type  2  customer  is  at  the  head  of  the 
queue. 
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if  Jl  + J2  <  S 


if  N  <  T  and  Jl  +  J2  =  S 


if  N 


if  Jl  +  J2  1  S 
and  Ql  =  Q2  =  0 


if  Jl  +  J2  =  S  and 
first  customer 
in  queue  is  Type  1 


if  Jl  +  J2  =  S  and 
first  customer 
in  queue  is  Type  2 


if  Jl  +  J2  1  S 
and  Ql  =  Q2  =  0 


if  Jl  +  J2  =  S  and 
first  customer 
in  queue  is  Type  1 


if  Jl  + J2  =  S  and 
first  customer 
in  queue  is  Type  2 


FIGURE    2.2 
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C.   STRUCTURE  OF  THE  ALGORITHM 

After  the  operation  of  the  system  was  understood  and 
represented  by  a  probability  model,  the  next  step  was  to 
write  out  a  simulation  algorithm  in  flow-chart  form.   The 
algorithm  must  be  capable  of  performing  at  least  the  follow- 
ing functions:   (i)  it  must  automatically  choose  an  appro- 
priate value  for  each  sojourn  time  based  on  its  current  proba- 
bility distribution,  (ii)  it  must  determine  which  of  the 
four  possible  events  is  to  occur  next,  (iii)  it  must  update 
the  state  of  the  system  based  on  the  type  of  event  that 
occurred,  the  value  of  Jl,  J2,  Ql,  Q2,  and  the  type  of  cus- 
tomer at  the  head  of  the  queue,  if  formed,  and  finally 
(iv)  it  must  gather  appropriate  statistics  on  the  figures 
of  merit.   The  general  flow-chart  of  Figure  2.3  describes 
the  way  in  which  the  program  performs  these  four  functions. 
Notice  that  the  algorithm  essentially  consists  of  one  loop 
embedded  in  another.   The  inner  loop  performs  the  major 
tasks  of  choosing  events  and  sojourn  times  and  updating  the 
system  variables.   This  loop  continues  to  iterate  for  a 
specified  length  of  simulation  time,  for  example  one  hour. 
When  the  time  limit  is  reached,  control  is  passed  to  the  outer 
DO  loop  which  gathers  interim  statistics  and  then  checks  to 
see  if  an  iteration  counter  has  reached  its  maximum  value. 
If  it  has  not,  then  control  is  passed  back  to  the  inner  loop 
and  the  process  starts  anew  with  the  values  of  the  system 
variables  preserved  from  the  end  of  the  previous  iteration. 
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GENERAL  FLOW  CHART  FOR  SIMULATION  ALGORITHM 
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If  the  iteration  counter  has  reached  its  upper  limit,  then 
final  statistics  are  gathered  and  the  algorithm  terminates. 
In  this  manner  one  day  of  system  operation  can  be  simulated 
by  setting  the  upper  limit  for  the  inner  loop  to  one  hour 
and  the  upper  limit  for  the  outer  loop  to  twenty-four.   The 
algorithm  uses  four  sub-algorithms  called  R0UTE1,  R0UTE2, 
R0UTE3  and  ROUTE 4 .   The  sub-algorithms  are  utilized  when  an 
event  occurs  and  their  function  is  to  correctly  update  the 
system  variables  according  to  the  rules  depicted  in  Figure  2.2 
Therefore,  for  instance,  R0UTE2,  which  is  utilized  when  a 
Type  2  arrival  occurs,  would  first  check  whether  all  the 
servers  were  busy;  if  not,  it  would  increase  the  number  of 
Type  2  customers  in  service  by  one.   If  the  servers  were  all 
busy,  R0UTE2  would  then  test  whether  the  number  in  the  system 
equaled  or  exceeded  the  queue  threshold;  if  not  then  the 
number  of  Type  2  customers  in  the  queue  would  be  increased 
by  one.   If  the  number  in  the  system  were  greater  than  the 
queue  threshold,  there  would  be  no  change  in  the  state  of 
the  system  and  control  would  be  passed  back  to  the  main 
algorithm.   The  flow-chart  for  R0UTE2  is  shown  in  Figure  2.4. 
The  flow-charts  for  the  other  sub-algorithms  are  similar. 

The  composition  of  the  queue  is  maintained  in  a  vector 
of  one's  and  two's,  called  NEXT1,  in  which  a  one  in  the  i 
position  indicates  that  the  i   customer  in  the  queue  is 
a  Type  1.   Each  time  a  customer  arrives,  another  element  is 
added  to  the  vector;  either  a  one  or  a  two  depending  on 
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the  type  of  customer  who  arrived.   Each  time  a  service  event 
occurs,  a  counter  is  incremented  by  one.   When  this  counter 
is  used  as  the  subscript  of  NEXT1,  it  points  to  the  element 
in  the  vector  which  represents  the  customer  currently  at 
the  head  of  the  queue.   For  example,  assume  that  five  arrival 
events  have  occurred  and  that  the  composition  of  the  vector 
NEXT1  is  as  shown  in  Figure  2.5.   Let  the  value  of  the 
service  counter  be  four  indicating  that  the  fourth  service 
event  is  about  to  occur.   Then  NEXT1(4)  indicates  that  the 
customer  at  the  head  of  the  queue  is  a  Type  2. 

Once  the  problem  was  understood  and  the  flow-chart 
written,  the  next  step  was  to  choose  an  appropriate  computer 
language  and  program  the  algorithm. 

D.   THE  COMPUTER  PROGRAM 

Fortran  IV  was  chosen  as  the  language  to  implement  the 
simulation  algorithm.   Initially,  consideration  was  given  to 
simulation-oriented  languages  such  as  SIMSCRIPT  and  GPSS, 
but  it  was  felt  that  their  sophistication  and  abstraction 
would  tend  to  obscure  some  of  the  inner  workings  of  the 
queuing  system  and  would  thus  offer  less  insight  into  the 
problem.   Coding  of  the  program  was  a  straightforward  trans- 
lation of  the  flow-chart.   A  program  listing,  with  explana- 
tory comments  and  a  sample  output,  is  provided  at  the  end 
of  the  thesis. 

The  three  figures  of  merit  are  printed  out  after  each 
hour.   These  are  cumulative  statistics  -  that  is,  they 
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reflect  the  operation  of  the  system  from  the  start  of  the 
simulation,  and  not  just  from  the  previous  hour.   The  printing 
out  of  cumulative  statistics  at  the  end  of  each  hour  was 
necessary  in  order  to  ensure  that  the  simulation  had  been 
allowed  to  run  long  enough  to  reach  the  practical  equivalent 
of  steady-state.   Note  that  there  is  little  fluctuation  in 
the  values  of  the  figures  of  merit  toward  the  end  of  the 
program  run.   If  there  had  continued  to  be  significant  fluc- 
tuation, it  would  have  indicated  that  the  system  was  not 
in  steady-state  and  that  the  program  should  be  run  again  for 
longer  simulation  time.   The  program  keeps  track  of  the 
composition  of  the  system,  by  customer  type,  at  all  times 
and,  for  each  customer  who  joins  the  queue,  it  computes  the 
actual  time  he  waits  for  service.   Thus,  the  average  waiting 
times  are  computed  by  simply  totaling  all  of  the  waiting 
times  and  dividing  by  the  number  of  customers  who  joined 
the  queue.   The  probability  that  there  are  no  customers  in 
the  system  is  actually  the  long-run  fraction  of  time  that 
the  system  was  in  that  state.   Similarly,  the  probability 
that  a  Type  2  customer  is  blocked  is  the  long-run  fraction 
of  time  that  the  queue  threshold  was  exceeded. 

The  program  utilizes  the  IMSL  sub-routine  RANDOM  to 
generate  the  pseudo-random  numbers  needed  for  determining 
the  sojourn  time  and  choosing  the  next  event.   RANDOM  takes 
as  input  a  "seed"  number  and,  for  a  given  seed,  will  always 
produce  the  same  string  of  numbers  which  are  uniformly 
distributed  on  the  range  from  zero  to  one. 
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Looking  at  the  print-out  from  the  simulation  program, 
the  first  question  that  might  be  asked  is:   "Are  the 
results  meaningful  —  do  they  accurately  represent  the 
functioning  of  the  real  system?"   Thus  the  need  for  some 
kind  of  verification  of  the  simulation  program  becomes 
obvious.   The  process  of  verification  is  detailed  in  the 
next  chapter. 
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Ill-   VERIFICATION  OF  THE  SIMULATION  ALGORITHM 

A.  NEED  FOR  VERIFICATION 

Verification  of  a  simulation  model  by  testing  it 
against  the  results  obtained  from  a  mathematical  model  is 
the  final  step  in  its  development.   As  mentioned  in  Chapter 
I,  the  verification  process  has  been  separated  from  the 
rest  of  the  development  of  the  algorithm  in  order  to  better 
highlight  some  interesting  results.   Of  course,  as  is  often 
the  case,  the  reason  for  going  to  a  simulation  approach  in 
the  first  place  was  that  no  general  mathematical  model 
existed.   However,  it  was  felt  that  if  a  mathematical  model 
could  be  developed  for  a  special  case  and  the  simulation 
verified  against  it,  then  it  would  at  least  indicate  that 
the  simulation  was  performing  properly.   This  is  the  approach 
that  was  taken,  and  it  is  detailed  below. 

B.  METHOD  OF  VERIFICATION 

1 .   The  Mathematical  Model 

It  was  decided  to  attempt  to  model  the  special  case 
of  the  problem  in  which  the  service  times  of  the  Type  1  and 
Type  2  customers  are  equal.  All  other  aspects  of  the  prob- 
lem remain  unchanged;  that  is,  the  arrivals  occur  according 
to  a  homogeneous  Poisson  process  with  different  parameters 
for  each  type  of  customer,  and  Type  2  customers  are  blocked 
from  joining  the  queue  when  the  queue  length  equals  or 


28 


exceeds  a  given  threshold.   Comparison  between  the  mathe- 
matical model  and  the  simulation  model  was  made  on  the  basis 
of  the  three  figures  of  merit  mentioned  in  Chapter  I;  those 
were,  (i)  the  probability  that  the  system  was  empty,  (ii)  the 
expected  waiting  time  for  any  customer  and  (iii)  the  probab- 
ility that  a  Type  2  customer  would  be  blocked.   The  equations 
for  (i)  and  (iii)  are,  although  complicated  and  tedious, 
straightforward  solutions  of  the  balance  equations  which 
resulted  from  application  of  the  theory  of  birth-death 
processes.    The  derivation  of  these  equations  is  detailed 
in  sections  two  and  four  of  this  chapter. 

It  appeared  that  the  simplest  way  of  computing  the 

2 
expected  waiting  time  was  to  utilize  Little's  equation  ; 


E[W]   -   — —   .  (3.d 


Where 


E [W]   =   the  expected  waiting  time  of  any 
customer 

E[Q]   =   the  expected  number  of  customers  in 
the  queue 

X  =   the  system  arrival  rate 


Gaver,  D.P.,  and  Thompson,  G.L.,  Programming  and 
Probability  Models  in  Operations  Research,  p.  467-482, 
Brooks/Cole,  1973. 

2Hillier,  F.S.,  and  Lieberman,  G.J.,  Introduction  to 
Operations  Research,  p.  291-292,  Holden-Day,  1967. 
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Little's  equation,  as  shown  above,  assumes  a  constant 
arrival  rate   and  is  not  directly  applicable  to  the  problem 
under  study  since,  as  was  discussed  in  Chapter  I,  the 
system  arrival  rate  depends  upon  the  number  of  customers 
in  the  system.   Therefore,  it  was  decided  to  modify 
Little's  equation  as  follows: 


Where : 


E[W]   =   SIS!  (3.2) 

E 


\        =   the  "effective"  rate  at  which  customers 
join  the  system. 


T-l 

X„   =  X-,    +  X-   y  p  (3.3) 

E       12  Ln    n 
n=0 


p   =   the  probability  that  there  are  exactly 
n  customers  in  the  system. 


Before  equation  (3.2)  can  be  verified,  it  will  be 
necessary  to  set  up  the  balance  equations  and  derive  an 
expression  for  p  .   At  the  same  time,  an  expression  for  the 
probability  that  the  system  is  empty  will  be  derived. 


Hillier  and  Lieberman,  p.  291-292 
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2.   A  Derivation  of  an  Expression  for  the  Probability 
That  the  System  is  Empty 

The  theory  of  birth-death  processes  provides 

expressions  for  long-run  or  steady-state  probability  that 

the  number  of  customers  in  the  system  is  a  given  value. 

To  determine  these  probabilities,  it  is  first  necessary 

to  define  some  terms.   Let 

N  =  the  number  of  customers  in  the  system, 

p  =  the  probability  that  N  =  n, 

X,  =  the  arrival  rate  for  Type  1  customers, 

X2  -  the  arrival  rate  for  Type  2  customers, 

y  =  the  service  rate  for  both  types  of  customers, 

a  =  overall  system  arrival  rate  when  N  =  i, 

T   =   queue  threshold;  that  is  when  N    T,  Type  2 
customers  are  blocked  from  the  queue, 

s   =   the  number  of  servers, 

y     =   system  service  rate  when  N  =  i. 


.  4 
It  can  then  be  shown  that  an  expression  for  p   is 

n 


.  (0)  .  (1)  .  (2)     (n-1) 

_    A *A »A .  .  .  a /  t  a  \ 

pn   "  "in    m    m       (hi     po  [J-q) 

y    •  y    •  y    •  •  •  y 


Gaver  and  Thompson,  pp.  467-482 
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which  comes  directly  from  the  basic  balance  equation. 
This  leads  to  the  following  expression  for  p  : 


P0   "   oo    i   (j-1}  (3.5) 

1  +  It   n  -n-r — } 
i=l  j=i^> 


The  remainder  of  this  section  is  devoted  to  manipulating 
equation  (3.4)  into  a  form  suitable  for  use  with  the  problem 
at  hand. 

The  first  step  is  to  write  down  the  system  arrival 

and  service  rates  that  result  from  the  queuing  discipline: 


J°>     »   +  i  (0)   -  n 

X  —       X ,  +  A  n  y      -   l) 


1    'v2 

1  +  X2 


X(1)   =   X,  +  X.  u(1)   =   u 


(2) 


1  (T_1)         5    +   1 

X  -   Al  +  X2 


(s) 

U       -   Sy 


x(T)   =   X 

x        Al  (s+1) 


y  Sy 


x(n»   =   x  P(n)   -  sp  . 
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Rewriting  equation  (3.5)  in  terms  of  the  above  expression 
yields  an  expression  for    in  terms  of  the  system  service  rate 
and  the  rate  at  which  customers  actually  join  the  system. 
This  is  an  important  point  and  it  will  be  repeatedly 
emphasized  in  the  remainder  the  the  thesis.   In  a  "normal" 
queuing  system,  the  arrival  rate  is  identical  to  the  rate 
at  which  customers  join  the  system.   However,  a  distinction 
must  be  made  in  any  system  in  which  the  queuing  discipline 
allows  customers  to  be  blocked  from  the  system.   The  correct 
expression  for  pQ,  in  terms  of  the  rate  at  which  customers 
join  the  system,  is  shown  below: 


P 


0  S  (\,+A9)n     T   (X.+A0)n         -     (X.+X9)T\,  (n  T' 

1+  [  y     l    2    +     y       1  ,2 — . —  +     y 

L,  .  n       L . ,   .   n-s   n     k 


(n-s)  n    ^i,   .  (n-s)  n 
n=l   n!y      n=s+l  s!s      y    n=T+l  s!s      y 


Next  examine  each  of  the  summations  in  the 
denominator  to  see  if  they  can  be  summed  in  closed  form. 

s  U+A2)n 

a.         J   ^—    •       This    expression   cannot   be    summed 

n=l      n!y 

in  closed  form;  however,  the  sum  can  easily  be  computed  in 
two  ways;  (i)  by  direct  summation  if  s  is  not  large  or 
(ii)  by  reference  to  a  table  of  cumulative  Poisson  proba- 
bilities.  The  second  method  is  possible  since 

((X  +X  )/y)n         A  +A2 
exp(-   ) 


n!  y 
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is  the  density  function,  with  respect  to  counting  measure, 

for  a  Poisson  distribution  with  parameter  — = -  .   Thus 

y 

if  the  figure  obtained  from  a  cumulative  Poisson  table 

look-up  with  the  appropriate  parameter  is  dived  by 

Xl  +  X2 
exp(-  )  then  the  desired  sum  will  result. 

T    (A1+X2)n 

b-         /_.. \ — —  .   This  sum  can  be  expressed 

n=s+l  s ! s      p 


in  closed  form  as  follows: 


— r  y         ,x    2     _      y         1    2    ] 

s!L  L  (n-s)  n      L  (n-s)  nJ 


n=s+l  s      u     n=T+l  s 


i,(Ws+1   :  .W.j 
5'-     Sus+1      jio     Sw 
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(x,+a0)         -   \,+x0  • 

(T+l-s)   (T+l)   £    l  Sy   M 
s         (j       j  —  u 


r'VV,s+1>     <w(T+1)    , 


s!  (1- 


A,+A0  u   a  (s+1)         (T+l-s)   (T+l) 

1   2,     Sy  S         y 

sy 


,.     ,T   (n-T) 
(A  +A  )  A 

c.         — - — t ; .   This  sum  can  also  be 

L ,  -,  ,     (n-s)  n 

n=T+l   sis      y 
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(xi+^2)T  h 1 


s!     x  (T+l-s)   (T+l) ;  v 

s        u        2. - 


Combining  the  above  sums  yields  the  desired 


expression  for  pn 


P 


°        a  U,+x 2)n    i     ,     a,+a9     U,+X,) (T+1) 

1+  V  x   2    |  r  x  (    -1    )  r  1   2    _    1   2 ,  . 

n=l   n!un      sl  \    h  +  X2      su(s+1)    s  (T+1"s)  „  (T+1) 


(3.6) 


^1+^2)T         Xl        1 

+  [ (    sl    } (  (T+l-s)  T+l} (   xT) 
s        u  -, 1 

Sy 


The  next  step  is  to  derive  equation  (3.2)  , 

the  modified  form  of  Little's  equation. 

3 .   Derivation  of  Little's  Equation  for  the  Case  Where 
the  Service  Rate  is  a  Function  of  the  Number  of 
Customers  in  the  System 

The  validation  of  equation  (3.2)  proceeded  in  two 

steps;  first,  the  equation  was  derived  in  general  to 

be  applicable  to  any  queing  system  and  next,  the  general 

equation  was  modified  to  yield  a  form  directly  applicable 


35 


to  the  problem  at  hand.   Before  beginning  the  derivation, 
it  was  necessary  to  define  a  number  of  quantities. 

T   =   the  queue  threshold. 

A (n)   =   combined  rate  at  which  customers  join  the 
system  when  N  =  n. 

A  =   the  event  that  a  customer,  of  either  type, 
arrives  and  joins  the  system. 

W   =   waiting  time  for  any  customer. 

p   =   probability  that  N  =  n. 

n  J 

In  the  general  case,  the  desired  form  of  the 
expression  for  the  expected  waiting  time  is: 

total  expected  waiting  time  for  all  customers 

_,  rwl      who  join  the  system  during  time  t 

total  expected  number  of  customers  who  join 
the  system  during  time  t. 


3.7 


In  order  to  write  the  above  expression  in  mathe- 
matical form,  first  consider  the  expected  waiting  time  of 
a  customer,  given  that  a  customer  arrives  and  joins  the 
system  when  N  =  n.   This  can  be  seen  to  equal  zero  if  there 
are  any  idle  servers  and  to  equal  the  expected  service 
times  of  all  customers  ahead  of  him  in  the  queue  plus  the 
expected  service  time  of  the  customer  being  served  if  the 
servers  are  all  busy  and  a  queue  is  formed.   This  can  be 
expressed  as: 
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E[W|A,N=n]       =     - 


if      n    <    s 


n-S  +  1  r; 

if      n    >    s 


Now  consider  the  expected  contribution  to  the  total 
expected  steady-state  waiting  time  which  results  from  a 
customer  arrival  in  a  short  time  interval  (dt) .   This  is 
simply  the  expression  derived  above,  multiplied  by  the 
steady-state  probability  that  a  customer  joins  the  system 
during  a  time  interval  dt  when  there  are  n  customers  in  the 
system  (A(n)dt),  and  by  the  probability  that  there  are  n 
customers  in  the  system  (p  ) ,  and  then  summed  over  the 
range  of  N.   The  expression  appears  as: 


1 

{^±i  A(n)p  }  dt   . 
L  su        n 

n=s 


If  this  last  expression  is  integrated  over  the 
range  (0,t)  then,  in  the  limit  as  t  goes  to  infinity,  the 
result  is  equal  to  the  total  expected  waiting  time  for  all 
customers  who  joined  the  system  during  time  t,  which  is  the 
numerator  of  equation  (3.7).   That  is: 


St  l  ^Hif±iA{n)Pn}dt   =     l   {2iri  A(n)pn}t  • 

0    n=s    y  n=s     M 
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The  denominator  of  equation  (3.7)  is  simply: 


I     (A(n)p  )t  . 
n=0       n 


Combining  the  expression  for  the  numerator  and 
denominator  yields  the  desired  expression  for  the  expected 
waiting  time  which  is  applicable  in  general: 


E[W] 


r  rn-s+l ,  ,  .    , 
L    H z Mn)p  } 

«— „   su       n 

n=s 

I    A(n)p 
n=0 


Equation  (3.8)  can  be  used  to  determine  the  expected 
waiting  time  at  steady-state  for  any  queuing  system  in 
which  the  rate  at  which  customers  join  the  system  is  a 
function  of  the  number  of  customers  in  the  system.   In  order 
to  make  equation  (3.8)  more  directly  applicable  to  the 
particular  system  under  study,  it  is  necessary  to  define 
some  more  terms.  They  are: 


co,  (n)   =   the  probability  that  a  Type  1  customer 
can  join  the  system  upon  arrival  when 
there  are  n  customers  in  the  system. 

co?  (n)   =   the  probability  that  a  Type  2  customer 
can  join  the  system  upon  arrival  when 
there  are  n  customers  in  the  system. 
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Note  again  the  important  distinction  made  between 
the  rate  at  which  customers  arrive  (A,  and  X_)  and  the 
rate  at  which  customers  joint  the  system  by  either  join 
directly  into  service  or  joining  the  queue  (A(n)).   The 
two  rates  are  related  as  follows: 


A(n)   =   x1u)i(n)  +  X  w2(n)  (3.9) 

Combining  equations  (3.8)  and  (3.9)  yields  the 
desired  expression  for  the  expected  waiting  time. 


CO  CO 

r    rn-s+l             ,       .  r    rn-s+l       ,    ,       , 

X.       I    {  — oj,  (n)p    }    +    X_       )     {— 0)o(n)p    } 

1      _        sy           In  2     z        sy        2          n 

E[w]       =      2_s nzs (3.10) 

00  CO 

Xl      ^n    wl(n)pn   +  X2      ln    w2(n)pn 

n=0  n=0 


Equation  (3.10)  can  now  be  applied  directly  to  the 
problem  of  interest  by  inserting  appropriate  values  for 
X,,Xp,uj,  (n)  and  oj_  (n)  .   However,  in  order  to  display  the 
connection  between  equation  (3.10)  and  equation  (3.2),  the 
modified  form  of  Little's  equation,  some  additional  work 
is  necessary. 

Recall  equation  (3.4)  in  section  two  of  this  chapter: 


,  (0)  .  (1)  .  (2)  .  (n-1) 

A      'A      •  A    •  •  •  A /  -3  a  \ 

Pn   "  II) TT) (T)        (nl    P0  *  (J*4) 

y    •  y    •  y     •  •  •  y 
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Writing  equation  (3.4)  in  recursive  form  yields; 


A 


Now  apply  the  definitions  above  to  substitute  for  X 


P 


S  "  Pn+1 


n      X  .  to,  (n)  + X  ..  w ~  ( n ) 


Now  insert  this  expression  for  p   into  the  numerator 


of  equation  (3.10); 


v  rn-s+l   /  x   i  ,    r  rn-s+l   ,  >   , 

X,     {  — to,  (n)p  }  +    \0      I    (  — co0(n)p  } 

1  L         Su    1     n  2      L         sy    2     n 

n=s  n=s 


n-s+1        S  y  Pn+1 


1        Sy     1     A,u,  (n)  +  X-to-ln) 
n=s  11        z  z 


n-s+1   ,_,    S  u  Pn+1 


+  X2        I     {~1T   U,2(n)X1a 

n=s  11        2  2 


(n-s+1)  ai1  (n)p   , 

y  r     ,   -1    n+1 
X,  uj.  (n)  +  X  oj  ( n 
n=s   1  1        2  2 


(n-s+1)  oo0  (n  p  ., 
V  r     2 n+1 

Xn  go,  (n)  +  X.i 
n=s   11        2 


X2  _A_  ^oi,  (n)  +  X0uo0(n)  ' 
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m 


Letting  m  =  n+1  yields  the  following  expression 
for  the  numerator  of  equation  (3.10): 


(m-s)  o)1  (m-1)  pm  oo        (m-s)oj  (m-l)p 

Al  mJ+1  X^dn-1)  +  X2u,2(m-1)  +  X2   Jg+1  x^dn-1)  +  ^/(m- 


»     A.ai.(m-l)  +  X  u)  (m-1) 
m=Ll   ^^^-D  +  X2.2(m-1)  (m"s)Pm 


I   (m-s)p  (3.12) 

m=s+l 


Notice  that  equation  (3.11)  is  an  expression  for 
the  expected  number  of  customers  in  the  queue  in  steady- 
state.   This  allows  us  to  rewrite  equation  (3.9)  as  follows 


E[W] 


E[Q] 


E[W 


X        I   Wl(n)p     +    x2      I   »2(n)p 

n=0  n=0 


_   expected  number  in  the  queue 

(expected  value  of  the  rate  at  which 
customers  join  the  system) 


Thus,  equation  (3.10)  can  be  seen  as  a  generalization 
of  Little's  equation  (3.2),  as  would  be  anticipated.   To 
complete  the  derivation  it  is  only  necessary  to  insert  values 
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for  w  (n)  and  a  (n)  as  appropriate  for  the  particular 

problem.   In  the  present  case,  these  values  are 


co,  (n)   =   1.0  for  all  values  of  n 


u)2  (n) 


—  i 


1.0   if   0  <  n  <  T 


0.0   if   n  >  T 


Inserting  these  values  into  equation  (3.9)  yields 

E[Q] 


E  [W] 


T-l 

1     2  Ln    n 
n=0 


which  is  identical  to  equation  (3.2). 

4 .   Calculations  of  the  Expected  Waiting  Time  and  the 
Probability  That  a  Type  2  Customer  is  Blocked 

Having  derived  equation  (3.2)  as  the  expression  for 

the  expected  waiting  time,  it  now  remains  to  obtain 

expressions  for  the  expected  number  in  the  queue  ( [EQ] )  and 

the  expected  rate  at  which  customers  actually  join  the 

T-l 
system  (X,  +  \  P  ) •   Once  again,  these  calculations 

n=0 
are  straightforward  but  tedious  and,  hence,  they  are  presented 

in  outline  form.   First  consider  the  following  expression 

for  the  expected  number  of  customers  in  the  queue: 
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E[Q]   =    I     (n-s)p 
n=s 


y  np    -  s    y  p 

„_„   n         n 
n=s        n=s 


Both  of  the  summations  above  can  be  expressed  in  closed 
form.   Details  of  the  necessary  calculations  are  included 
in  Appendix  A  and  only  the  results  of  those  calculations  are 
given  below: 


x1+x2 

-tQ]   =  Sft^'l-^tV-^   "7^   ] 

( 1 — = — -)   ( 1 — - — -)  ^ 

sy  sy 


U,+A„) 


(T+l)  Al+A2 


[         T+l     +       sy    j  j 


(T+l-s)(T+l)  L    A.+A0        An+A   9 

sy         sy 


s!     s(T+l-S)u(T+l)1   V      V2 

sy        sy 


r     SPq     r/l+A2,s    (A1+X2)(T+1) 

1      A,+A0  U   y   '      (T+l-s)  (T+l) 
1   2,  s       y 


s:(1-fr^ 


(A1+,\-)T         A. 


+  *Po'   Si       ,T+l-sl  (T+l)  -\"-        (3-13' 

s       y      -, _ 

sy 
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Next,  consider  the  expression  for  the  expected  rate 
at  which  customers  actually  join  the  system; 


T-l 

A    =   A,  +  A    }   p 
E       12  L-    pn 
n=0 


Al  +  V?o  +   ?  Pn] 


n=l 


*,  +  *0tP  +  p  I  — - — — 

1     2  ^o   ^o    ,    ,  n 

n=l   n!U 


T-l  (A  +  A  )n 
P°  n=s+l  s-s(n-s)un 


Both  of  the  summations  above  were  encountered  in 
deriving  the  expression  for  p_;  the  first  sum  cannot  be 
expressed  in  closed  form,  but  the  second  can  and  details  of 
the  calculations  can  be  found  in  section  two  of  this  chapter 
The  final  form  for  A   is : 


s  (A,+A,Jn 


n=l   n! u 


.       (A+A  )(S+1)    (A  +A  )T 

+  r ± ( ± t — ± ± — ]}       (3.14) 

A  +A   (    (s+1)  (T-s)  TJJ 

s.(l-^_£)   s^  ^ 
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The  final  figure  of  merit  is  the  probability  that 
a  Type  2  customer  is  blocked  from  joining  the  queue  (P~  ) . 
This  is  simply  the  probability  that  the  number  in  the 
system  equals  or  exceeds  the  queue  threshold  which  is: 


P2b   =    L  Pn      1  P~  ' 


T-l 

D 

n 


n=T  n=0 


T-l 
Since  the  summation,  I      p   was  just  encountered  in  finding 

n=0   n 
A£  it  is  possible  to  write  the  desired  expression  for  P  , 

directly : 


s  (VVn 

P-.        =      1   -   p    {1   +      I   — -± S 

2b  *o  L.       ,    n 

n=l   n!y 


i       (W(s+1)    <vv 


+  [ nrx-( b+tj (t^t^]}  (3-15) 

al(l-   1     2)         s    y  s  y 

sy 


5.   Results  of  the  Validation  Efforts 

At  this  point,  expressions  have  been  derived  for  all 
three  of  the  figures  of  merit;  the  probability  that  the  system 
is  empty  (Pn)/  the  expected  waiting  time  for  any  customer 
(E[W]),  and  the  probability  that  a  Type  2  customer  will  be 
blocked  from  joining  the  system  (?2h^  *   T^e  next  steP  i-s  to 
compute  values  for  these  quantities  and  compare  the  results 
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obtained  with  the  values  of  those  same  quantities  which 
result  from  running  the  computer  simulation. 

In  order  to  reduce  the  variance  in  the  simulation 
results,  the  method  of  antithetic  variables  was  used.5 
Recall  from  Chapter  Two  that  the  RANDOM  subroutine  utilized 
by  the  program  produces  a  string  of  pseudo-random  numbers 
which  range  from  zero  to  one  and  that  it  will  produce  the 
same  string  each  time  unless  the  "seed"  number  is  changed. 
In  the  antithetic  variable  method  of  variance  reduction  two 
simulation  runs  are  made  for  each  data  point  desired.   In 
the  first  run  the  random  numbers  are  used,  unmodified,  to 
determine  the  sojourn  time  and  choose  the  next  event.   In  the 
second  run,  the  program  is  modified  slightly  so  that  the 
random  numbers  are  subtracted  from  one  forming  an  antithetic 
random  variable  which  is  then  used  as  before.   Results  from 
the  two  runs  are  then  averaged  and  the  value  becomes  the 
data  point  for  that  iteration.   The  seed  number  is  then 
changed  to  produce  a  new  random  number  string  and  the  process 
is  repeated  for  as  many  iterations  as  desired.   Finally,  the 
data  points  for  all  of  the  iterations  are  averaged  to  produce 
an  estimate  of  the  variable  under  consideration.   As  illus- 
trated in  Gaver  and  Thompson  [Ref .  4] ,  sampling  by  antithetic 
variables  can  reduce  the  variance  in  simulation  results  by 


Gaver  and  Thompson,  pp.  584-58  6. 
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more  than  one-half  of  the  variance  expected  using  straight- 
forward sampling  techniques. 

Finally,  the  comparison  of  simulation' and  analytical 
values,  which  serves  as  the  validation  of  the  simulation  model, 
is  presented  in  Figure  3.1.   The  values  shown  under  the 
simulation  heading  were  obtained  by  averaging  over  the  results 
of  ten  iterations  of  the  antithetic  variable  sampling  tech- 
nique just  described;  that  is,  twenty  simulation  runs  were 
made  in  all.   Each  run  simulated  thirty-six  hours  of  system 
operation.   The  comparison,  measured  in  terms  of  the  percent 
difference  between  the  analytical  and  simulation  values, 
seems  to  offer  a  good  indication  that  the  simulation  model 
is  functioning  as  intended.   Notice  that  the  largest  discrepancy 
is  consistently  in  the  expected  waiting  time  values.   No 
explanation  for  this  evidently  systematic  effect  was  immediately 
apparent. 
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Verification  of  the  Simulation 


Case  I   A   =  1.0   A2  =  .5 


U1  =  1.0 


M2  =  1.0   p  =  .695 


Analytical 

E[W]      .58234 
pQ        .174387 
p2b       .220708 


Simulation 

.5721195 

.173633 

.219648 


Difference 

1.75 
.43 
.48 


Case  II  A   =  1.0 


A2  =  .9 


y1  =  i.o 


u2  =  i.o 


p  =  .799 


E[W] 


'2b 


.74334624 
.10332954 
.33665 


.757377 
.104199 
.3381565 


1.89 
.841 
.447 


Case  III 


X±   =    1.0 


A2  =  .99 


u1   =  i.o 


u2  =  i.o 


p  =  .816 


E[W] 


2b 


.77443362 
.09074268 
.36099 


.7874915 

.093319 

.3624095 


1.69 
2.839 
.393 


Case  IV  A   =  1.75   A2  =  .249 


v1   =  i.o 


y2  =  i.o 


p  =  .913 


E[W] 


2b 


3.475469 
.043550164 
.69541 


3.58214 
.044149 
.695702 


3.07 
1.37 
.042 


48 


IV.   AN  ANALYTICAL  APPROXIMATION 

A.   RATIONALE  FOR  DERIVING  AN  ANALYTICAL  APPROXIMATION 

In  the  previous  three  chapters  a  simulation  model  was 
developed  and  verified.   It  was  seen  that  the  simulation 
appeared  to  do  a  reasonably  good  job  of  modeling  the  system; 
good  enough,  at  least,  for  most  engineering  purposes.   Why, 
then,  go  to  the  trouble  of  developing  a  second  model?   There 
are  several  reasons.   First,  all  of  the  reasons  for  developing 
the  simulation  apply  equally  to  an  analytical  model.   In 
particular,  the  process  of  developing  an  analytical  approxi- 
mation allows  an  opportunity  to  gain  additional  insight  into 
the  problem  from  a  different  perspective.   Such  insight  could 
lead  further  toward  the  ultimate  goal,  which  is  an  exact 
but  simple  and  useable  analytical  solution  if  it  exists. 
Secondly,  the  implementation  of  the  approximation  described 
in  this  chapter  requires  tools  no  more  sophisticated  than  a 
pocket  calculator.   Thus,  the  analytical  approach  offers 
some  practical,  financial  advantages  over  the  simulation. 
Financial  advantages  are  even  more  apparent  when  the  need  to 
do  sensitivity  analysis  is  considered.   In  order  to  determine 
the  sensitivity  of  the  system  to  a  change  in  one  variable, 
it  is  necessary  to  do  a  large  number  of  simulation  runs  for 
each  different  value  of  the  variable.   Consider  Figure  3.4 
in  which  it  was  necessary  to  do  eighty  simulation  runs  in 
order  to  be  able  to  compare  four  variable  changes.   The 
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analytical  approximation  need  only  be  done  once  for  each 
variable  change  and  if  done  on  a  computer,  the  calculations 
required  take  only  a  fraction  of  the  time  used  by  one  simula- 
tion run.   On  the  other  hand,  it  would  be  easier  to  modify 
the  simulation  to  investigate  other  figures  of  merit  than 
it  would  be  to  write  new  equations  for  the  approximation 
model. 

B.   DEVELOPMENT  OF  THE  ANALYTICAL  APPROXIMATION 

Recall  that  in  Chapter  Three  it  was  shown  that  an 
analytical  solution  to  the  problem  existed  in  the  special 
case  in  which  the  service  rates  were  equal  and  that  expressions 
were  developed  for  the  expected  waiting  time,  the  probability 
that  the  system  is  empty,  and  the  probability  that  a  Type  2 
customer  will  be  blocked.   This  suggested  approximating  the 
unequal  service  rate  case  by  taking  a  weighted  average  of 
the  service  rates  and  using  the  resulting  value  in  the  equa- 
tions for  the  equal  service  rate  case.   The  intended  weighting 
was 


u   =   Ul  hi +  y2        T~l (4-1} 

a.,  +\„   y   p  i  +  -,    i  p 

1   2  n-0   n       1      2   nio   n 


Where : 


y  =  the  weighted  average  of  service  rates 
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All  other  variables  are  as  defined  previously.   The  weighting 
factors,  which  are  the  ratio  of  each  type  of  customer's 
arrival  rate  to  the  overall  system  arrival  rate,  are  intended 
to  represent  the  long-run  proportion  of  time  that  each 
operator  would  spend  serving  each  type  of  customer.   Unfor- 
tunately, this  weighting  scheme  is  infeasible  since  it 
requires  values  for  the  p  's  which  are  not  available.   There- 
fore, an  iterative  procedure  was  decided  upon  in  which  another 
type  of  weighting  was  used  to  get  approximate  values  for  the 
Pn's  which  could  then  be  used  in  equation  4.1  to  get  better 
approximations  which  are  again  used  in  equation  4.1  and  so 
on.   The  initial  weighting  utilized  was 

The  rationale  behind  this  weighting  scheme  is  the  same  as 

before,  however,  it  can  be  seen  that  this  weighting  does 

not  take  into  account  the  effect  of  blocking  Type  2  customers. 

A  computer  program  was  written  which  would  perform  the 
necessary  calculations  and  iterate  the  procedure  as  many 
times  as  desired.   A  program  listing  and  sample  output  are 
provided  at  the  end  of  the  thesis. 

C.   VERIFICATION  OF  THE  ANALYTICAL  APPROXIMATION 

As  before,  the  final  step  in  developing  a  model  is  to 
test  its  performance  by  comparison  with  a  suitable  mathe- 
matical model.   Recall  that  the  simulation  model  was  verified 
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by  testing  it  against  an  analytical  solution  to  the  special 
case  of  the  problem  in  which  the  two  types  of  customers  have 
equal  service  rates.   This  would  suggest  trying  to  find  another 
special  case  of  the  problem  in  which  the  service  rates  are 
different  and  for  which  an  exact  solution  could  be  found 
with  which  to  verify  the  analytical  approximation.   Unfor- 
tunately, such  a  special  case  was  not  discovered.   An 
attempt  was  made  to  solve  the  problem  for  the  case  of  a 
single  server  with  the  queue  threshold  equal  to  one  where 
Type  2  customers  could  only  access  the  system  when  it  was 
empty.   Although  it  seemed  that  a  solution  should  be  obtain- 
able, the  complications  associated  with  the  different  service 
rates  made  even  this  comparatively  simple  problem  unsolvable 
in  the  time  available. 

Since  the  approximation  model  cannot  be  verified  analytically, 
the  only  alternative  available  is  to  verify  it  by  comparison 
to  the  simulation  model.   Of  course,  the  simulation  is  itself 
an  approximation  since  it  is  only  exact  in  the  limit  as  time 
goes  to  infinity  and  true  steady-state  is  reached.   However, 
it  was  shown  in  Chapter  Three  that,  at  least  for  small  systems 
and  equal  service  rates,  the  results  obtained  from  running  the 
program  for  thirty-six  hours  of  simulation  time  differ  from 
the  exact  analytical  solution  by  less  than  three  percent. 
Chapter  Two  described  how  the  simulation  program  was  able 
to  keep  track  of  the  identity  of  the  customer  at  the  head 
of  the  queue  and  thereby  correctly  adjust  the  system  service 
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rate  when  a  state  change  occurred.   Therefore,  the  simula- 
tion should  be  equally  valid  in  the  case  where  the  service 
rates  are  different.   Thus,  a  reasonable  amount  of  agree- 
ment between  the  approximation  and  the  simulation,  while  not 
conclusive,  should  be  sufficient  to  provide  confidence  in 
the  approximation  model. 

To  see  how  closely  the  two  models  agree  in  two  test 
cases,  refer  to  Figure  4.1.   The  same  method  of  antithetic 
variable  sampling  that  was  discussed  in  Chapter  Three  was 
used  to  obtain  the  results  for  the  simulation  model.   The 
approximation  results  were  obtained  by  carrying  out  ten 
iterations  of  the  procedures  discussed  earlier  in  this  chap- 
ter.  It  is  interesting  to  note  that  the  largest  change  in 

the  figures  of  merit  occurred  between  the  first  iteration  — 

— ?  ^1        ^2 

which  utilized  \i       =  y,  •? — 77 —  +  y~  -r — —j—   —  and  the  second 

iteration,  and  that  there  was  little  change  between  succeeding 

iterations.   In  fact,  there  was  no  change  at  all  after  the 

fifth  iteration.   This  effect  is  apparent  in  the  sample 

output. 

The  degree  of  agreement  between  the  two  models  shown  in 

Figure  4.1  indicates  that  either  model  would  be  appropriate 

for  a  first  estimate  of  system  performance.   The  choice  of 

model  would  therefore  depend  on  the  relative  merits  of  each 

as  evaluated  in  a  particular  situation.   For  instance,  if 

sensitivity  analysis  was  to  be  performed,  the  approximation 

model  would  be  appropriate,  while  if  an  estimate  of  the 


53 


Analytical  Percent  Difference 

Approximation  Simulation         ( (Sim- Anal) / (Sim) 


Case  I    X      =  1.0  X2    =  .5   u1  =  .666   u2  =  .833 


s  =  3   T  =  6 


E[W]  .45110463  .4700145  4.02 

p  .1047012  .099507  5.22 

p_  .10048  .1064845  5.63 
^2b 


Case  II   \.  =  1.75   A2  =  .249  \i±   =  1.0   y2  =  .833 

s  =  3   T  =  6 

E[W]     .39178724            .40591  3.47 

.10854438            .109263  .657 

.11799               .12258  3.74 


P0 


P2b 


Figure  4.1 
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proportion  of  customers  whose  waiting  time  was  greater  than 
a  given  value  was  desired  then  the  simulation  model  should 
be  chosen. 
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V.   CONCLUSION 

Two  aspects  of  this  problem  which  made  it  interesting 
were  the  different  service  rates  for  the  two  types  of 
customers  and  the  fact  that  Type  2  customers  are  blocked 
from  the  system  whenever  the  queue  threshold  is  reached. 
Recall  that  Chapters  One  and  Two  discussed  the  implications 
of  these  aspects  of  the  problem  and  showed  that  they  made 
an  exact  analytical  solution  difficult,  if  not  impossible, 
to  achieve.   In  particular  it  was  shown  that  the  different 
service  rates  required  that  any  solution  be  able  to  keep 
track  of  both  the  number  and  type  of  customers  in  the  queue 
and  in  service  as  well  as  the  identity  of  the  customer 
at  the  head  of  the  queue.   Also,  the  blocking  of  Type  2 
customers  implies  that  the  systems  arrival  rate  is  not  con- 
stant but  is  a  function  of  the  number  of  customers  in  the 
system. 

The  only  real  attempt  made  to  find  an  analytical  solu- 
tion utilized  the  theory  of  birth  and  death  processes  and 
it  is  possible  that  some  other  analytical  approach  could 
yield  results.   However,  given  the  complications  mentioned 
above,  it  seems  that  any  future  effort  expended  on  this 
problem  could  best  be  directed  toward  improving  the  approxi- 
mate solutions. 

It  would  be  worthwhile  to  attempt  to  verify  both  the 
simulation  and  the  analytical  approximation  with  an  exact 
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solution  to  a  special  case  with  unequal  service  times. 
Chapter  Four  mentioned  that  such  a  verification  was  attempted 
but  abandoned  primarily  due  to  time  constraints.   It  is 
quite  possible  that  additional  work  in  this  area  could  be 
successful. 

The  simulation  algorithms  should  be  studied  to  determine 
confidence  limits  on  its  estimates  and  to  determine  the 
sensitivity  of  the  algorithm  to  changes  in  system  parameters. 
This  type  of  analysis  would  also  serve  to  increase  the  con- 
fidence in  the  comparison  between  the  simulation  and 
approximation  models. 

Finally,  there  is  obviously  much  room  for  improvement 
of  the  analytical  approximation.   Although  the  model  developed 
in  Chapter  Four  appears  to  perform  fairly  well  and  has  an 
intuitive  rationale,  there  exist  a  number  of  techniques 
which  could  be  expected  to  yield  better  results.   In  particu- 
lar the  theory  of  Gaussian  approximations  and  diffusion 
approximations  discussed  by  Gaver  and  Lehoczky  [Ref.  2  and  3] 
and  Kleinrock  [Ref.  6]  should  work  well  at  least  in  the 
case  where  the  system  is  heavily  saturated. 
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APPENDIX   A 

DETAILS    OF    THE    DERIVATION    OF    THE    EXPRESSION 
FOR    THE    EXPECTED    NUMBER    IN    THE    QUEUE     (E[QJ) 


E[Q]       =         I    (n-s)p        =         I      npn   -    s      £      Pr 
n=s  n=s  n=s 


«  T    (n)(X1+A2)np0  -       (n)(X1+X2)TX1n"TpQ 

lm  nPn      =         L   ,     (n-s)    n        +        *      ,    9n-s)    n 

n=s  n=s    s!sv         ;y  n=T+l        s's  'y 


T    (n)(X1+A2)np0  Pq      -    (n)(A1+X2)n 

a* }  (n-s)    n  s!  L„£„      m  (n-s)    n 


n=s      s ! s  y  n=s      s 


n(A.+X9)n 

I 


(n-s)    n 
n-r+1    s  y 


Pq      »    (n)(X1+A2)n  PQ^2S 

1* )       s!  (n-s)    n  s!  (       y       ] 

n=s   s  y 


2 

x,+x9  (x,+x2) 

s+(s+l)    -^7T  +   (s+2)  2    ,         + 

M  s   y 


(s+k)  (X   +X2)k 

+  iTk +    '••    ] 

s    y 
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£0      ?    (n)(VX2)n 
si  (n-s)    n 


n=s      s 


y 


_      P0fAl+A2.sf     ,         1         .    ,      "    A+A2,n 


_W    '  n=]/      sy  X_A1+ 


sy 


sy 


P0(V^2)S[ 


si  y 


+ 


Xl+X2 

sy 


1    Al+A2  (1.A1+A2 


sy 


sy 


pn        o°      (n)  (X-+X-) 

si         ,L,  (n-s)    n 

n=T+l      s  y 


n 


p0   (x1+x2) 


T+l 


s!  T+l-s    T+l 

s  y 


A-|+A,j                            A  ..  +  A  _     ~ 
[(T+l)  +  (T+2)    -= -+   (T+3)  (-= -)       + 


sy 


sy 


Xi+X2   k-l 

+     (T+k)   (-i -)K    L    + 

sy 


.] 


p0  (x1+x2) 


T+l 


s!       T+l-s    T+l 
s  y 


[  (T+l 


+ 


h+X2 

su 


A    +A      '           A    +A 
1 —  ...         (1 — = -) 


sy 


sy 


n 


?    (">(W    P0         P0r/l  +  A2,2r s_ 

(n-s)    n        '     s! 1 (      y      '     L      A,+A 


VX2 

sy 


n=s    s!s 


y 


A    +X 
1-    ...    ^       (1-    t.       ) 


sy 


sy 


A       A       T+l 
k     1       2}  r       T+l 


Al+A2 

sy 


T+l-s    T+l         A.+A_ 
s  y  -       1      2 

sy 


Xl+X2    2 

(l-  ;„  V 


sy 
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0°   (n)  (A1+A9)TA.n"Tpn 

b.)    J    )     2.     1 12. 

n=T+l   s!s(n-s)un 


s! 


nA 


n-T 


n=T+l  s(n"s)yn 


P0(A1+A2) 


T+l-s  T+l[(T+1)+(T+2)i^ 


Xl    2 

+  (T+l)  (-±P  +  ...] 


-   (n)(A1+A2)TA1n-TpQ 
n=T+l   s-s(n-S)yn 


P0(A1+A2) 
s! 


Xl  1 

± r  (T+l) — = — 

T+l-s  T+1L  K         '       A, 
s      M  , 1 

sy 


+ 


1  2 

(1 — -) 

sy 


I    np 


n=s 


n 


p0   Al+A2  s 


l- 


Xl+X2 

sy 


+ 


Al+X2 

sy 


Xl+X2  2 

(1-    i.    ) 


sy 


(A1+A2) 


T+l 


[   T+1, 


Al+X2 
sy 


T+l-s  T+lu   A.+A0 
s      y     ,  _  1   2 

sy 


A1+A2N2 


]} 


1- 


sy 


P0(A1+A2) 


s! 


Al  r       T+l 

T+l-s  T+l L.  Xl 

s     y    1 

sy 


^1 

sy 


Xl    2 
sy 
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2.       s      I   p 


n=s 


n 


.£|iW»o 


n=  s    sis 


(n-s) n 


,,     .,     >T,     n-T 
(At+A0)     Xn  p 


1    "2'      vl 


n=T+l    s!s(n-S)yn 


(A.l) 


n 


T    (A1+X2)"p0 
a .  )       s      I   — 


,     (n-s)    n 
n=s    sis  y 


spn     (\+X,  )S      -      X    +X  (A    +  X    )T+1      oo      X    +X 

L  s  .  ^    sy       ;  T+l-s    T+l     >n         sy    '     J 

y  3=0  s  y  3=0 


s! 


s   p 


0 


x  +x 

s! (1-    i,    Z) 


[( 


Al+A2,s 


(x1+x2) 


T+l 


T+l-s    T+l 
s  y 


sy 


b.)       s        I 


,,     ,  .     x  T,    n-T 

(X  +x  )    xx        p0 


n=T+l   s-s(n-s)yn 


S?0[ b! ( 


T+l-s    T+l 

s  y 


sy 


I    P 


s   p 


0 


n=s 


n 


X   +X 

s!  (1-    _,,    Z) 


Al+A2    s 


(X  +x2) 


T+l 


T+l-s    T+l 

s  y 


sy 


+  po[ el ( 


xl  1 
± )  ( — ± — ) 

T+l-s    T+l;  v      xn ' 
s  y  , 1 

sy 


(A. 2) 
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Finally,  combining  A.l  and  A. 2  yields  the  result 
shown  in  Chapter  3. 
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FRCGRAM  LISTING  FQo  SIMULATION  ALGORITHM 


C 

C 

c 

c 
c 
c 

c 

r 

c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


THIS    PROGRAM    SIMULATES    A    MULTI-SERVER    QUEUING    SYSTEM 
WHICH    IS    BEING    ACCESSED    BY    TWO   TYPES    OF    CUSTOMERS.    THE    TWO 
TYPES    OF    CUSTOMERS    HAVE    CIFFERENT    ARRIVAL    AND    SERVICE 
RATES.    TYPE    1    CUSTOMERS    ARE    ALLCWED    TC    JCIN    THE    S\STEM    RE- 
GARDLESS   OF    QUEUE    LENGTH.    TYPE    2    CUSTOMERS    ARE    BLOCKED 
FROM    THE    QUEUE    WHEN   THE     QUEUE    LENGTH    EXCcECS    A    GIVEN    ThRES 
HOLD.    THE    SIMULATION    IS     EASED 
TIMES    BETWEEN    EVENTS    WITH    THE 
CHOSEN    RANOCMLY. 


ON    RANDCNLV    CHOSEN    SOJOURN 
NEXT    EVENT    TO    OCCUR    ALSC 


12 
15 


20 


27 


10 


THE  VARIA8 
LAMDA1 
LAMDA2 
MU1=SE 
MU2=SE 
KQ=THE 
KS=THS 
ol  =  THE 
J  2- THE 
NUM1=T 
NLM2=T 

NBLCCK=THE 
I=THE 
J  =  THE 
SOJTI M 


LES    USED    I 

=ARRIVAL    F 

=ARFIVAL    R 

RVICE    RATE 

RVICE    RATE 

QUEUE    THF 

NUMBER    OF 

NUMBER    OF 

NUMBER    OF 

HE    NUMBER 

HE    NUMBER 

NUMBER    OF 
NUMeER     CF 
NUMEER    OF 
=THE    SCJCL 
1/(LAM0A1 


N    THE    PROGRAM    APE    DEFINED    EELOVs. 
ATE    OF    TYPE    1    CUSTOMERS 
ATE    OF    TYPE    2    CUSTOMERS 

OF  TYPE  1  CUSTOMERS 

OF  TYPE  2  CUSTOMERS 
ESHOLD  FOR  8LCCK  ING  OP  TYPE  2 

SERVERS 

TYPE    1    CUSTOMERS    IN    SERVICE 

T>P?    2    CUSTOMERS    IN    SERVICE 
CF    TYPE    1    ARRIVALS 
OF    TYPE    2    ARRIVALS 

TY°E    2    CUSTCM5FS    BLOCKED. 
BOTH    TYPES    WHC    HAVE    JCINSC    THE    QUE 
BOTH    TYPES    WHO    HAVE    LEFT    THE    QUEUE 
FN    TIME=AN    EXPONENTIAL    RANDOM    VARI 
♦LAMDA2+J1*MU1+J2*MU2) 


PE 

CO 

1UV 

CI 

1ME 

RE 

FO 

RE 

FO 

WR 

FO 

WP 

FO 

1  'A 

ISC 

1  'T 

l'C 

1,T 

WR 

FO 

l'C 

1'2 

l'Q 

16/ 

CA 

SA 

IF 

Jl 

J2 

NU 

NU 

cu 

KU 
KU 
TI 
2E 

ei 
i  = 
j= 


AL*4  LAMCA1 
MNCN  LAMDA1 
C1,NUMQ2,NU 
MENS  I  CM  U(2 
(2000) 

AC (5, 5)  LAM 
RMAT(5F10.5 
AC(5,1D  )  Jl 
RMAT(7I5  ) 
115(6,12)  K 
RMAT(T45, » N 
ITE(6,15J  L 
RNAT(T45,' I 
RRIVAL  RATE 
,F10.5, /,T2 
>FE  2  SERVI 
N  time*  ,T50 
2 C,1  QUEUE  T 
ITE(6,20J  J 
RMATU45,'  I 
F  TYPE  1  IN 
IN  SERVICE 
LSLE'  ,T50,I 
/) 

UU  OVFLCW 
VE2=LAM0A2 
(ITEM. EG. 1) 
=  C 
=  0 

MC1  =  0 
MC2=0 
MV»  1=0.0 
MNQD=0 
MIN=0 
ME  =  0.  0 
PTIM=0.0 
KTIM=0.C 
0 
C 


,LAMDA2  ,MU1,MU2 
,L  A^IDA2t  TIME,  STIME  ,QT 
M1,MUM2,MBL0CK, I, J »N0 
),NEXT1 (2000)  ,QTIMS(2 

CA1,LAMCA2,MU1,MU2,TM 
) 

,J2,NUMQlTNUMQ2f  KSTOF 


STOP 
UMBER    0 
AM  DAI, L 
NITI  AL 
•  ,T50,F 
C,'TYPE 
CE   RATE 
,F10.5, 
FRESHCL 
1,J2,NU 
MTIAL 
SERVIC 
•,T50,I 
8,/,T20 


F  ITERA 
AM0A2,M 
CONCITI 
10.5,  /, 
1  SERV 
',T50,F 
/,T20,' 
C' ,T50, 
MQltNUM 
STATE  V 
E«  ,T50, 
8,/,T20 


TIONS  T 
U1.MU2, 
ONS'  ,// 
T20, »TY 
ICE  RAT 
10.5,/, 
NUMBER 
18//  ) 
Q2 

ARIABLE 
I3,/,T2 
,'  NUMBE 
R    CF    TY 


IMS ,KQ,KS,I >,J1,J2,N 

kAIT  ,NEXT1 

CCC ),STIME(20CO) ,UTI 

A> 

,KS,KQ 


H  IS    RUN'  ,14//) 
1MAX,KS,KQ 
,T20,»TYFE    1    ', 

PR    2    ARRIVAL    RATEST 
E'  ,T50,F1C5,  /,T20, 
T20,  'MAX.SI^UATI'  , 
CF    SERVERS'  ,T50,I8  ,/ 


S*  ,//,T20, ' NUMBER  ' 


R  Qc 
FE  2 


NUMBER  OF  TYPE 
TYPE  1  IN  ', 
IN  QUEUE'  ,T5C, 


GO  TC  27 


DO  500 
NUM1=0 


K=1,KST0P 
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NUN2=C 

NBLCCK=0 

NOWAIT=0 

NLINC=NLMG1+NUMQ2 

30  K2-NUMQH-NUMQ2 
IF(K2.GS.KQ)     GC    TC    3000 
LAMCA2=SAVE2 

CO    TO    3001 
3CCC       LANC*2=0.0 

3001  D=LAMDA1+  LA MD A2+ J  1* MU 1+ J  2*MU2 

C  THE    FCUR     FOSSIBLE    EVENTS    ARE    TYPE    1    ARRIVAL?     TYPE    2 

C    ARRIVAL, TYPE    1    SERVICE    AND    TYPE    2    SERVICE.    THE    EVENTS    ARE 
C    CHOSEN    RANDOMLY    AND    THE    NECESSARY    CHANGES     IN    THE     STATE 
C    VARIABLES    Jl, J2tNUMGl,NUMQ2    ARE    MADE    BY    CALLING    ThE    AFFRC- 
C    PRIATE    SUBROUTINE. 

CALL    RAND  CM  (IX,U,3) 

S0JTIM  =  (-1.0/C)*ALCG(U  (1)  ) 

IF(K2.LT.KC)     GO    TO    2002 

BLKTIM=BLKTI M+SOJTIM 

T2INT=(-1.0/SAVE2I*AL0G(U(3  )) 

IF( T2INT. LE  .SCJTIM)     NB LCCK=NBLOCK 4 1 

3002  TIME=TIME+SOJTIM 
*SYS=J1+J2+NUNC1+NLNC2 

IF  (KSYS.GT..1  )    GO    TC    3003 
IE  FTIM=ZERTIM  +  SOJTIM 
C 

c 

3003  IF(L( 2) .LE.LANDA1/D)    GO    TO    31 
IF(L(2).LS. (LAMDA1+LAMCA2)/C)     GO    TC    32 
IF<U(2>.  L£.  (LAMDA1+LAMD42+J1*MU1)/D)    GO    TO    33 

IF  (J2.EQ.0)    GC   TO    41 
CALL    RCUTE4 
GO    TC    49 

31  CALL    RCUTE1 
GC    TC    49 

32  CALL    R0UTE2 
GO    TO    49 

33     IF  (Jl.EO.O)     GC    TO    41 

CALL    R0UTE3 

GO    TC    49 
4i    WRITS<6,42) 
42    F0RVAT(T55  ,  'ERROR. SERVICE    ATTEMPTEC    kITH    0    CUSTOMERS1) 

STCF 

49  FK=FLOAT(K) 
IF(TIM£.LE.FK*TMAX)  GO  TO  30 
t*RITE(6,50)  KtTIME 

50  FCFNATI//T45, 'END  OF  INTERV AL », 13 ,/, T45 ,  • CUMULATIVE  ', 
I'TINE  IS'  ,F12.5) 

C 

C 

C  THE    FOLLOWING    BLOCK    CF    STATEMENTS     (l)COMPUTES    Tl-E    mjn- 

C    EER    UFO    QUEUED    (2JC0MPUTES    THE    WAITING    TIME    FOR    E4CF 

C     (WTIVE(M)    ANC     (3IC0MPUTES    THE    AVERAGE    MITING    TINE     (CUN- 

C    AV  G  )  . 

NCC=I-NLINQ 

NLlNC=NUMQi+NLNQ2 

WTSLM=0.0 

DO    8C    M  =  l  ,J 

VTI^E<M)  =  STIME(M)-GTIME  (M) 

WTSLN=WTSUM+WTIME (M) 
80    CONTINUE 

CUKWT=CUMWT+WTSUM 

KUMNQD  =  KUVNCD-»I 

KUMIN=KUMIN+NLM1+NUM2 

CUMAVG=CUMWT/<KUMIN-(I  -J)  » 

WRITE(6»90)     CUMAVC- 
90  FORfAT(/T20t 'CUMULATIVE    AVERAGE    WAITING    TIME     FCP    ALL     ■ 

1, •CLSTOMERS* ,/,T20,  'EXCEPT    THCSE    BLOCKED1  ,  F12  .  I ) 

FC=ZERTIM/TIME 

URIT£(6,92)    PO 


64 


92 


93 


2C0 
250 


F0RNAT(//T20t  'CUMULATIVE    PROPCRTICN    CF    TIME     THAT    N  =  0  ' 
1F12.5) 
PKC  =  BLKTIM/TIMS 
WRITE<6,92)     PKC 

FQRNAT(//T20,  'CUMULATIVE    PRQPGP.TICN    CF    TIME    THAT    TYPE 

BE    BLCCKED'  ,F12  .5  ) 


500 


1,'     I    CUSTOMERS' 
IF(  (I-J) .EQ.O ) 
ISTCF= I-J 
CO    200    L-l.ISTOP 

CTIfE(L  )=CTIME(J  +  L) 

CONTINUE 

1  =  1- J 

J  =  0 

CONTINUE 

STOF 

ENC 


t/  ,T2C,1 WOULO 
GO    TO    250 


SUBROUTINE    LISTINGS    FCR    SIMULATION    ALGORHI-N 
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SUBROUTINE    RCUTE1 

THIS    SUBROUTINE    IS    CALLED    IN    EVENT    CF    A    TYPE    1    ARRIVAL 

REAl*4    LAMDA1,LAM0A2,MU1,MU2 

COMNCN  LAMCA1 , LAMDA2 ,TI ME , STI ME , QT INS tKO ,KS,IX,J1,J2,N 
1LMQ1,NUMQ2,NUM,NUM2,NEL0CK,I  ,J  ,NCW4  17,NEXT1 

CI  MENS  I  ON  U(2) t NEXT  1 ( 2 CCO ) » QTIME ( 2CC D ) ,STIME( 2300 ) ,WTI 
iw=  (2C00) 

NUN1=NUM1 +1 

IF  (J1  +  J2.EQ.KS)     GO    TO    10 

Jl=Jl+l 

NOWMT^NOfcAIT-U 

GO    TO    20 

NUNG1=NUMQ1+1 

1  =  1  +  1 

NEXTKI  )=1 

CT^EI  I)=TIME 

RETURN 

END 


2C 


30 


SUB 
THI 
RE  A 
COM 

1UMQ 
DIM 

1NS( 
NUM 
IF  ( 
J2  = 
NOW 
GO 
NUN 
1  =  1 
NEX 
QTI 
FET 
ENC 


ROUTINE    RCUTE2 

5    SLBROUTINE    IS    CALLED    IN 

L*4    LAfDAl,LAMDA2»ML'l,MU2 


EVENT    CF    A    TYPE    2    ARRIVAL 


NCN    LAM0A1,LAMDA2,TIME,STIM€,QTIME,KQ ,KS  t  IX, Jl, J2t N 

1,NUMQ2,NUM1,NUM2,NBL0CK » I , J,NCWA I T, NEXT1 

EN SI  ON    U(2) ,NEXT1 (2  000), QTIME (2000)  »ST  IME( 2000), WTI 

*>  /*  ^  «  ■ 


^lCO) 

2=NUM2+1 

J1+J2.EQ.KS) 

J2  +  1 

AIT=N0V»AIT+1 

TO    30 

C2  =  NLMC2-H 

+  1 

IK  I)  =  2 

*E<  I)  =  TIME 

LRN 


GO    TO    20 


SUBROUTINE    R0UTE3 

THIS    SUBROUTINE    IS     CALLED    IN    EVENT    OF    A    TYPE     1    SERVICE 
RcAL*4    LAMDA1,LAMDA2,MU1,MU2 

CCMNCN    LAMDA1,LAM0A2,TIME,STI NE , GTI NE , KQ,  KS, I  X , Jl t J  2 ,N 
1LMQ1  ,NUMQ2  ,  NU M , NUM2 , NBLOCK ,1  ,  J  ,NCHA  IT  ,  NE  XT1 
CI  MENS  I  ON    U(2),NEXT1(20C0),QTIME(2CCC)  ,  ST  IME  (  20CO  )  ,WTI 
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1MS(  2000) 

KCUE=NUMQ 1+NUMQ2 

IF( KQUE.GT..1J    GO 

J 1  =  J  1-1 

GO    TC    30 

1C 

J  =  J+1 

1F(  NEXTKJ)  .GT.l) 

MJPC1=NUMQ1-1 

GO     10    25 

20 

Jl=vl~l 

J2=J2+1 

NUVG2=NUMQ2-1 

25 

STIME (J) =  TIME 

20 

RETIRN 

ENO 

TO    10 


GO   TO    20 


10 


20 

25 
30 


SU8FCUTINE    R01TE4 

THIS    SUBRCUTIME    IS    CALLED    IN    EVENT    CF    A    TYPE     2    SERVICE 

FE4l*4    LAMCA1 ,LAMDA2 »PU1» MU2 

COMMON  LAVCA1?LAMDA2,TIMEtSTIMP,QTINEtKQtKS,I>t  J1,J2,N 
1UMC1»NLMQ2»NUV1 ,NUM2 ,NB LOCK , I , J , NCW k  IT ,NEXT  1 

CIMENSICN  U(2  JtNEXTl  (2000) ,QT I NE ( 2CC C  ) , STI ME  (  20 CO )  , fcT I 
1VE(2000) 

KQU  =  =  MIMQ1  +  NUVC2 

Ic(KQUE.GT..l )    GO    TC    10 

v2  =  o  2-1 

GO    TC    20 

J=J  +  1 

IF(NSXT1(  J). GT.l )    GC    TO    20 

MIVC1=NUMQ1-1 

J1  =  J1+1 

J2=J2-1 

GO    TC    25 

NjyC2=NUMQ2-l 

STI  fS(  J)  =TIME 

RETLRN 

ENC 
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PROGRAM    LISTING    FOR    THE    ANALYTICAL    APPROXIMATION 

C  THIS    DOUBLE    PRECISION    FRCGRAM    COMPLIES    (1)    A    WEIGHTED 

C   AVERAGE     SYSTEM    SERVICE    RATE    (DMU),(2)THE    PROBILI7Y    THAT 
C    THE    SYSTEM    IS    EMPTY    (P0),(3)THE    EXPECTEC    NUMBER    CF    CUSTO- 
C    NERS    IN    THE    QUEUE    (EQ),AND    (4)THE    EXPECTED    WAITING    TIME 
C    FOR     4NY    CLSTOMER     ( EW  )  .    THE    GENERAL    EQUATIONS    OEVELCFEO    IN 
C    CHAPTER    2    OF    THE    THESIS    ARE    USED    IN    AN    ITERATIVE    PROCEDURE 
C    AS    DESCRIBED    IN    CHAPTER    4.    TO    IMPROVE    READABILITY,     THE    EQ- 
C    LATICNS    ARE    BROKEN    DOWN    INTO    A    NUMBER    OF    TERMS    (TERM    THRU 
C    TERM8) . 


20 
25 


100 


IMPLICIT 

DIMENSION 

1EW  (20)tP2 

RE*C(2,13 

13  FQRMT<40 

kR  ITE(6,5 

5C  FQRMAT(4D 

FEAC  (2t 18 

18  FORNATOI 

WRI 7E(6t5 

55  FOPNATOI 

CO    1000    I 

IF(I  .GT.l 

CMU (I )=(J 

GO    TC    25 

CMU( I) =(U 

SUM=O.OD 

DO    10 

SUM1  = 

CCNTI 

Fl*1.00O- 

F2=(DL1+D 

F3=3.0CO- 

TE  PM  1=1.0 

TERN2=(  (C 

1(( C11+DL2 

TERM3=( (D 

1CML <I)**( 

P0( I  1*1.0 

TEP*4=< (C 

1(( DL1+DL2 

PSUMI) =p 

CLEFFd  )  = 

TERM5=(  (  D 

TEF^6=(  (C 

KKT  +  1)  )  )* 

TEPM7=  ((P 

KKTU-KS  ) 

1 (KS*CML(I 

TE  PN8=(  (  ( 

1(KS**(KT  + 

EQ( 3)=(P0 

1(TERM8-»TE 

EW(I)=EQ( 

P2B(I  )=1. 

WPITE  (6,2 

200         FORNAT(// 

1C15.8  ,T10 

1000       CONTINUE 

STCF 

ENC 


REAL*8     (A-HtO-Z) 

DMU(2C)  , PC(2C) ,PSLM(20)  , CLEFFC20  )  ,£C< 20 ) 
B(  20  ) 

)DL1  ,DL2»L1,U2 
15.8) 

C)0LltDL2,LltU2 
15.8) 

)    KStKT,  IMAX 
6) 

5)KS  ,KT t  IMAX 
6) 

=  1,  IMAX 
)    GC    TO    20 
1*DL1+U2*0L2)/  (DL1  +  DL2) 


1*QL 

0 

0    J  = 

SUM1 

NUE 

(CL1 

L2  )/ 

CL1/ 

C0/( 

Ll+D 

)**( 

Ll  +  O 

KT  +  1 

D0/( 

Ll  +  D 

)**K 

0(1) 

DL1  + 

Ll  +  D 

Ll  +  C 

((KT 

0(  I) 

*0VU 

)> )/ 

CL1+ 
1-KS 
(I)/ 
RM3) 
I)/D 
ODO- 
00)0 
T5,D 
6,01 


1+U2*DL2*PSUM(I-1) )/<CLl+DL2*PSUMI-l) ) 

1    KS 

+  ! CL1  +  CL2)**J/(FACT(  J  MCMUd  >**J) 


+DL2 )/(KS 
<KS*DMU( I 
(KS*DMU(I 
FACT(  KS  )* 
12  )**(KS  + 
KT+1)  /(  KS 
L2 ) **KT*D 
)*(1.0C0- 
l.ODO+SUM 
12)**(KS+ 
T/  (KS**(K 
*(1.0DC+S 
CL2*PSUN  ( 
L2)/0MU(I 
L2)**( KT+ 

♦  D/F1+F2 

*(DL1+CL2 
(I )**  (KT  + 
F3**2  ) 
CL2)  /DMU( 
)*OMU  (I)* 
FACT ( KS)  ) 


*D*U( 
)  ) 
)  ) 
F  1) 
1  )/  (K 
**(  KT 
Ll)/( 
CLl/( 
1+TcR. 
1)/(K 
T-KS) 
LM1  +  T 
I  ) 

)  )  **K 
1  )/(K 
/Fl** 
)**KT 
D)  )* 

I)  )** 

*  (KT  + 

*(TER 


I)  ) 


S  * CMU  (  I ) 
+1-KS)*D 
FACT( KS) 
KS*DMU(I 
M1*TSRM2 
S*CMU(  I) 
*DMU( I  )* 
ERM4/( FA 

S*  (KS/F1 
S**<KT+1 
2) 

*DL1)/  (F 
((KT+1  )/ 

KS-(DL1+ 
1) ) )/(  FA 

M5-TERM6 


LEFFd  ) 

FSUM(I) 

MU(  I  )  tPOd  )tPSUM(I  )  ,  CL 

15.8,T26,D15.3  ,T46,D15 

2.5) 


**<KS+1)  )  )- 
ML(I)**(KT+1) ) ) 
*KS**(KT+1-KS)* 
)  )  )  ) 
HTERM3  ) 
**(KS+1) ) )- 
*KT)  ) 
CT(KS)*F1 ) ) 

+F2/F1**2 ) 
-KS)*DMU (I )** 

*CT(KS)*KS** 
F3+(DL1/ 

CL2  )**(KT+1)/ 
CT(KS)*F1  ) 
l+TERM7-KS*P0<  I)* 


EFF(I)  ,EM  1)  ,P2B(I  ) 
.8,T66,D15  .8,T86  , 


COUeLE  PRECISICN  FUNCTION  FACT(N) 
C      THIS  SLBROUTIN'E  USES  A  GAMMA  FUNCTICN 
C   FACTORIAL  OF  THE  INTEGER  N. 

IMPLICIT  RE4L*8(A-H,0-Z  ) 


TO  COMPUTE  THE 
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IF  (N.GT.l)    GO    TO    10 
FACIEI  .ODO 
RETURN 
10  K=N+1 

CN  =  CFLQAT (K  ) 

FAC1=CGAMMA(DN) 

RETURN 

ENC 
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S-ANFL5  OUTPUT  FOR  ANALYTICAL.  AP  PR  CX  I  f  ATION 
SYSTEN  FAFAMETERS 


TYPE  1  ARRIVAL  RATE  0.175000COC  01 
TYPE  2  ARRIVAL  RATE  0.249CCCCCD  GO 
TYPE  1  SERVICE  RATE  0.10CCOOCOD  31 
TYPE  2  SERVICE  RATE  0.833C00C0D  00 
NUMBER  OF  SERVERS  3 
CUEUE  THRSSHOLC  6 
NUMBER  OF  ITERATIONS  10 

ITERATICN  NUMBER   1 

THE  WEIGHED  AVG.  SERVICE  RATE  IS 
THE  V.4LUE  OF  PO  IS  0.107874C9D  00 
THE  EFFECTIVE  SYSTEM  ARRIVAL  RATS 
THE  EXPECTED  WAITING  TIME  IS 
THE  FROB.  THAT  TYPE  TWOS  ARE 


0.97919£10D  00 


IS  0.1S6931220  01 
0.39563682C  00 
BLOCKED  IS  C.11922827C 


ITERATICN  NUMBER 


THE  W5IGI-TE0  AVG.  SERVICE  RATE  IS 
THE  VALUE  OF  PO  IS  0.10S55151D  CO 
THE  EFFECTIVE  SYSTEM  ARRIVAL  RATE 
THE  EXPECTED  WAITING  TIME  IS 
THE  FROB.  THAT  TYPE  TWOS  ARE 

ITERATICN  NUMBER   2 


0.9814C2C7D  00 


01 


ITERATICN  NUMBER   4 

THE  WEIGHTED  AVG.  SERVICE  RATE  IS 

THE  VALUE  OF  PO  IS  0.10854433D  00 

THE  EFFECTIVE  SYSTEM  ARRIVAL  RATE 

THE  EXPECTEC  WAITING  TIME  IS 

THE  FROB.  THAT  TYPE  TWOS  ARE 


0.98127888D  00 


01 


ITERATICN  NUMBER 


CC 


IS  0.19696232D 
0.39174657C  CO 
BLOCKED  IS  C.11797899C  00 


THE  WEIGHTEC  AVG.  SERVICE  RATE  IS  0.98127663D  00 

THE  VALLE    OF    PO    IS    0.1095443CD    00 

THE  EFFECTIVE  SYSTEM  ARRIVAL  RATS  IS  0.19696199D  01 

THE  EXPECTEC  WAITING  TIM£  IS  0.39178767C  00 

THE  PROB.  THAT  TYPE  TWOS  APE  BLOCKED  IS  C.11799220C  00 


IS  0.19696200D 
0.39178722C  00 
BLOCKED  IS  C.117992C6C  CO 


THE  WEIGHTED  AVG.  SERVICE  RATE  IS  0.98127888D  00 

THE  VALUE  OF  PC  IS  0.10854438D  00 

THE  EFFECTIVE  SYSTEM  ARRIVAL  RATE  IS  0.1969620Q0  CI 

THE  EXPECTED  WAITING  TIME  IS  0.39178724C  CO 

THE  PROB.  THAT  TYPE  TWOS  ARE  BLOCKED  IS  C.117992C6C  CO 

ITERATICN  NUMBER   6 

THE  WEIGHTED  AVG.  SERVICE  RATE  IS  0.96127888C  00 

THE  VALUE    OF    PO    IS    0.10854438D    00 

THE  EFFECTIVE  SYSTEM  ARRIVAL  RATE  IS  0.1S6962OOD  01 

THE  EXPECTED  WAITING  TIME  IS  0.39176724C  00 

THE  PROB.  THAT  TYPE  TWOS  ARC  BLOCKED  IS  C.117992C6C  CO 
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SANFLE  CUTPUT  FOR  SIMULATION  ALGORITHM 


NUMBER  OF  ITERATIONS  THIS  RUN   36 


TYPE  1  APRIVAL  RATE 
TYPR  2  ARRIVAL  RATS 
TYPE  1  SERVICE  RATS 
TYPE  2  SERVICE  RATE 
MAX. SIMULATION  TIME 
NUMBER  CF  SERVERS 
QUEUE  THRESHOLD 


INITIAL  CCNOITICNS 

1  .75000 
0.24SCC 
1.00000 
0.833GC 
60.00000 

3 

3 


NUMBER  CF  TYPE  1 

NUMBER  CF  TYPE  2 

NUMBER  CF  TYPE  1 

NUMBER  CF  TYPE  2 


INITIAL  STATE  VARIABLES 


IN  SERVICE 

IN  SERVICE 

IN  QUEUE 

IN  QUEUE 


END  OF  INTERVAL   1 

CUMULATIVE  TIME  IS     62.C3427 

CUMULATIVE  AVERAGE  WAITING  TIME  FOR  ALL  CUSTOMERS 
EXCEPT  THOSE  BLOCKED     -7.11838 


CUMULATIVE  PROPORTION  OF  TIME  THAT  N=0 


0.13202 


CUMULATIVE  PROPCRTICN  OF  TIME  THAT  TYPE  2  CUSTOMERS 
WOULC  eE  BLOCKED      0.55560 

END  OF  INTERVAL   2 

CUMULATIVE  TIME  IS    122.51188 

CUMULATIVE  AVERAGE  WAITING  TIMS  FOR  ALL  CUSTOMERS 
EXCEPT  TFCSE  3L0CKED    -14.22775 


CUMULATIVE  PROPORTION  OF  TIME  THAT  N=0 


0.06789 


CUMULATIVE  PROPORTION  OF  TIME  THAT  TYPE  2  CUSTOMERS 
WOULC  BE  BLOCKED      0.77146 
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END  OF  INTERVAL  33 

CUMULATIVE    TIME     IS       1982.23901 

CUMULATIVE    AVERAGE    WAITING    ttmE    FOR    ALL    CUSTOMERS 
EXCEPT    ThCSE    BLOCKEC      -234. 92163 

CUMULATIVE    PPOPORTICN    OF    TIME    THAT    N=0  0.00420 

CUMULATIVE  PPOPCRTION  OF  TIME  THAT  TYPE  2  CUSTOMERS 
WCULD  BE  ELOCKEC      0.98587 

END  OF  INTERVAL  34 

CUMULATIVE  TINE  IS   2041.77637 

CUMULATIVE  AVEFAGE  WAITING  TIME  FOR  ALL  CLSTOMERS 
EXCEPT  THOSE  BLOCKED   -242.04044 

CUMULATIVE  PROFCRTICN  OF  TIME  THAT  N=0      0.00408 

CUMULATIVE  PROPORTION  OF  TIME  THAT  TYPE  2  CUSTOMERS 
WOULD  BE  BLOCKED      0.98628 

END  OF  INTERVAL  35 

CUMULATIVE  TIME  IS   2101.31372 

CUMULATIVE  AVERAGE  WAITING  TIME  FOR  ALL  CUSTOMERS 
EXCEPT  THOSE  BLOCKED   -249.15926 

CUMULATIVE  PPOPCRTION  OF  TIME  THAT  N=0      0.00396 

CUMULATIVE  PROPORTION  OF  TIME  THAT  TYPE  2  CUSTOMERS 
WOULD  BE  BLOCKED      0.98667 

END  OF  INTERVAL  36 

CUMULATIVE  TIME  IS   2160.65107 

CUMULATIVE  AVERAGE  WAITING  TIME  FOR  ALL  CLSTOMERS 
EXCEPT  THOSE  BLOCKS   -256.27808 

CUMULATIVE  PROPCRTION  OF  TIME  THAT  N=0      0.00385 

CUMULATIVE  PROPORTION  OF  TIME  THAT  TYPE  2  CUSTOMERS 
WCULC  BE  ELOCKEC      0.98704 
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